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1 . Introduction 


Most of the problems in fluid mechanics are multi- 
dimensional and time-dependent as well. There are few excep)- 
tions, which can be expressed using only one time-like variable 
and one space-like variable. The most important are: spherical 
and cylindrical time-dependent flows, time-dependent flows in 
ducts of constant or slightly varying cross-sectional area, and 
steady, supersonic, two-dimensional or axisymmetric flows. We 
will generically call them "one-dimensional flows". In this pa- 
per, we will focus our attention on those problems where the 
time-like variable is really time. 

What are the reasons for limiting the present analysis to 
one-dimensional flows? What are the advantages and disadvan- 
tages? What is extendable to multi-dimensional flows and vdiat is 
not? We will begin by outlining some answers to these questions; 
such answers will guide us in the present study. 

The mathematics of one-dimensional problems is much simpler 
than the mathematics of multi-dimensional problems. From a 
numerical standpoint, this is reflected in limited storage, sim- 
ple coding, and the possibility of using a large number of compu- 
tational nodes. Consequently, numerical tests for accuracy can 
be conducted by doubling the number of nodes, over and over 
again. Discontinuities of any kind are perpendicular to the 
flow, so that their nature is as simple as possible. Finally, 
there is a certain number of realistic problems for which the ex- 
act solution is known. They can be used as benchmarks for a stu- 
dy of accuracy in an absolute fashion. Finally, basic ideas 
about procedures for accelerating a computation and devices such 
as multiple grid or multigrid techniques can be also studied. 

The main disadvantages of one-dimensional analysis proceed 
from the unrestricted and uncritical use of the above advantages. 
The mathematical simplicity of one-dimensional flows, indeed, re- 
flects the physical difference in complexity between one- and 
multi-dimensional flows. In the latter, waves propagate in more 
than one direction. One should welcome progress in one- 
dimensional problems with a guarded optimism towards a better un- 
derstanding of multi-dimensional problems. Extensions from one- 
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dimensional ideas and techniques to their multi-dimensional coun- 
terparts is far from being a matter of routine. Some methods 
which have been specifically created for one-dimensional flows 
may prove to be very difficult to extend to multi-dimensional 
flows. Methods which require a large number of nodes and whose 
stepsize is severely limited by the fineness of the mesh may 
become prohibitive in more than one dimension. Finally, discon- 
tinuities in multi-dimensional problems require a knowledge of 
the orientation of their normal, regardless of the computational 
technique being used. Nevertheless, jumps across such discon- 
tinuities (which obviously occur along the normal to them) can be 
investigated in a one-dimensional way. 


2. Historical background 

The goal of numerical calculations of unsteady flows is an 
accurate description of evolutions of such flows in time (use of 
unsteady codes to generate steady solutions asymptotically, a 
device which has been very helpful in the first two decades of 
electronic computing, is tapering off, for economical reasons). 
Accuracy implies a detailed description of wave propagation, in- 
cluding various degrees of discontinuities. A proper description 
of wave propagation relies on defining the domain of dependence 
of the wave and coding the discretized equations of motion accor- 
dingly, The importance of this idea was recognized (but, unfor- 
tunately, not emphasized) by Courant, Isaacson and Rees in an 
early work [1]. Godunov [2] went even farther, by suggesting to 
consider the local discretization of the Euler equations as the 
solution of a series of local Riemann problans. For many years, 
Godunov's suggestion lacked to find developers; some applications 
[ 3 ] oversimplify the original idea, and the accuracy of their 
two-dimensional applications seems not to be very high. 

More recently, integration techniques have been proposed, 
and are still under development, which emphasize the role of the 
domain of dependence and which show a proper concern for the 
analysis of discontinuities. The explicit schemes of Steger [4], 
Roe [5], van Leer [ 6 ] and Osher [7] as well as the X-scheme C 8 ], 
all belong to this category. Their importance stems from the 
possibility of describing complicated multi-dimensional flows. 
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with transonic effects and an evolving pattern of interacting 
discontinuities. 

The schemes of Steger , Roe, van Leer and Osher all derive 
from the Euler equations in divergence form. For each one of the 
equations, the pertinent flux is decomposed in parts which are 
carried by different waves. The waves are defined by the eigen- 
vectors of the system. For each partial flux, the proper domain 
of dependence is easily found and the discretization of space 
derivatives is carried on accordingly. All schemes can be made 
to have second-order accuracy. Special tests must be performed, 
and modifications must be made to the standard schemes, to assure 
monotonicity. Shocks are captured and smeared out over three in- 
tervals; contact and gradient discontinuities are also captured 
and smeared out over an increasing number of intervals as the 
computation proceeds (we make these statements on the basis of 
results presented by Pandolfi [9]. who uses a variant of Osher’ s 
scheme) . 

The x-scheme proceeds from the Euler equations written in 
terms of Riemann parameters and entropy. The definition of 
domains of dependence is obvious, the coding is extremely simple, 
and second-order accuracy is easy to achieve. Shocks appears as 
junps over a single interval, but they may not move at the right 
speed or lock in the right position, and they do not provide any 
jump in entropy, unless special provisions are taken. Contact 
discontinuities are captured properly and smeared out as in 
Osher ’s scheme. Monotonicity at strong gradient discontinuities 
is not assured. 


3. An improved x-scheme 

We intend now to show that very accurate solutions of un- 
steady, one-dimensional compressible flow problems can be ob- 
tained using a coding technique inspired by the x-scheme, plus 
adequate tracking of discontinuities. We will show that the 
handling of discontinuities can be coded, once and for all, 
essentially using Boolean algebra. The management of the discon- 
tinuous features of the flow is, thus, of no concern for the 
user. In this respect, the computational technique presented in 
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this paper is equivalent to the other techniques mentioned above 
[4, 5, 6, 7, 9 ]. We consider our technique more efficient for the 

following reasons: 

1) All discontinuities are captured within one single 

mesh . 

2) The basic integration scheme is much simpler; conse- 
quently, it is also faster. 

3) The additional logic, related to the presence of 
discontinuities, is treated in the fastest possible way, through 
the use of Boolean algebra. 


4. Basic equations 

We will denote velocity, density, pressure, speed of sound 
and temperature by u, p, p, a, and T, respectively. In a non- 
dimensional form, T=p/p, a^=yT. In addition, define m=pu, 
e=p/(Y-1 )+mu/2, h=u(p+e), Prln p, and s=(P-y In p)/[y(y-D] 

(non-dimensional entropy) . 


The one-dimensional Euler equations in divergence form are 

= 0 ( 1 ) 

where 
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- 
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p+mu 
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( 2 ) 


Alternatively, the equations can be written in the form 
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In characteristic form, (3) are written as follows: 


2P^ + X,(p 


r*^x1 a x1 


u„) t X2(P,2 = 0 


-i 


“xl' * ^2*’’x2 * a “x2* ' ° 


=t * u = 0 


where 


= u-a, X^ = u+a 


( 4 ) 


(5) 


( 6 ) 


As explained in C8], the subscripts x1, x2 and x all denote x- 
derivatives. In discretizing (5), however, the derivatives 
denoted by xi (i=1,2> are approximated by finite differences tak- 
en from the side from which X^ proceeds. The derivative denoted 
by X alone is approximated by a finite difference taken from the 
side from which u proceeds. 


Riemann invariants can easily be introduced in (5) if the 
flow is isentropic. In this case, indeed, 

dP = i i® (7) 

6 a 


where 6 = (y-1)/2 for brevity. Therefore, letting 
the first two equations (5) beccxne: 

I ®t + ^1^1x1 ^2^2x2 = ° 

^ ‘^t ” ^1^1x1 ^2^2x2 ■ ° 


( 8 ) 


(9) 
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whereas the third equation (5) becomes an identity. If the flow 
is not isentropic, (7) must be replaced by 


dP.x|a-,ds 


( 10 ) 


and (9) by 

i ®t " 2^®t ^2^^2x2 ■ ®®x2^ " ° 

^ “t - - ^=x1> * '‘2<"2X2 - = ° <"> 

. u Sx = 0 

Note that, to maintain the spirit of (5), in the first two equa- 
tions (11) the finite differences approximating the derivatives 
of s are taken in the same direction as the finite differences 
approximating the derivatives of that to which they are con- 
nected. In fact, it is not the domain of dependence of a or s 
separately which counts, but the domain of dependence of P. The 
derivative, s^ , in the first equation (11) is obtained from the 
third equation, where s^ is approximated as specified for (5). 

If the flow in a duct of variable cross-section, A=A(x), is 
considered in a quasi-one-dimensional way, (11) can still be 
used, provided that a term, 2auA^/A is added to the left-hand 
side of the first equation. 


5. Basic integration scheme 

The basic integration scheme for (11) is a variation of the 
X-scheme proposed independently by Zhu et al. [10], slightly 
modified by Pandolfi and Colasurdo (Politecnico di Torino, priva- 
te cOTimunication) and reformulated by the senior author of this 
paper in a manner which seems easier to extend to multi- 
dimensional problems [11]. 
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Let 


® ®x1^ 


^2 - - ^2^^2x2 “ ^ ®x2^ 


^3 = - “ "x 


= - a u Ajj/A 


Eqs. (11) can thus be written in the form: 

= I (f^ + f2 + 2af3 + 2f^) 


“t = 2 ^^2 “ ^1^ 


( 12 ) 


( 13 ) 


s 


t 


f 


3 


Let n denote the point to be computed (x=nAx) and k the time in- 
dex (t=kat). Assuming, for argument's sake, that u is positive, 
f 2 , f^ and ^4 are approximated by 


^2= - 


^ 2 .n~^ ^ 2 ,n -1 *^ 2 .n~ ^ 2 ,n- 1 ~ ^n^^n~^n- 1 ^ 

2 AX 


( 14 ) 


^3 * " 2 


“n^ “n -1 V ®n -1 


AX 


^4 “ - ^n^n^V^^n 


and f^ is approximated by 

^l.n'*' ^ 1 ,n +1 ^ 1 ,n+ 1 ~ ^ 1 .n~ ^n^^n+ 1 ~ ^n^ 


"l 


if X^< 0 , and by 


AX 


^Kn"*" ^ 1 ,n -1 ^ 1 .n~ ^ 1 ,n- 1 ~ ^n^^n" ^n- 1 ^ 
^1 “ ■ 2 AX 


( 15 ) 

(16) 


( 17 ) 


( 18 ) 


if X^> 0 . 
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To simplify the coding, we define X^=u and R-=s. The in- 
tegration step is divided into a predictor and a corrector level. 
At the predictor level, we go through four steps. Each step is 
executed at all computational nodes: 

1) Computing and (8), and X^ from (6), and 

letting x^ =u; computing the stepsize according to the CFL condi- 
tion. 


2) Computing R.^ (i=1,2,3) 

way that R. =(R. ._r. )/Ay 

IX, n i,n+1 i,n'^' ' 

3) Defining 


and storing them 


in such a 


j=n - ^ tsgnCX^^^) +1] (19) 

to determine the domain of dependence of R. and using (14), (15) 
and (17) after replacing n with j, and (16) as is; denoting any 
fj^ so obtained by f^, and storing them in a double array. 




i,n 




4) Using f^ to evaluate (13) and then 

X fH 


k+1/2 k k .. 
®n = «n ^ «t,n"^ 


( 20 ) 


( 21 ) 


for g = a,u,s. Note that (21) performs an update at At/2 only, 
because of the factor 1/2 in (20). 


Boundary conditions, when needed, are enforced between step 
3 and step 4. They are discussed in Section 13. 


At the corrector level, similar 
ceeds without changes till (20), but 


steps are taken. One pro- 
(20) itself is replaced by 


^k+ 1 /2 -D -k 
i,n ■ 1 ■ ^i,m 


( 22 ) 


where 


m=n-sgn(X^^^) (23) 

and (21) remains formally the same in the coding but in practice 
will mean 
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k»1 , k,1/2 ^ gk*l/2jt 
n n t.n 

A straightforward but tedious analysis shows that the pro- 
cedure outlined above has second— order accuracy even with non- 
constant coefficients [12]. Stability, diffusion and dispersion 
are easily studied in the isentropic case. By addition and sub- 
traction of the two equations (9), we obtain: 


^1t ^1^1x1 = ° 

^2t * ^2®2x2 ° 


(24) 


The two unknowns are now uncoupled. The two equations express 
the well-known fact that is conveyed along ( u-a)-characteris- 
tics and R 2 along ( u+a)-characteristics. The analysis can be 
carried through for each one of the equations, separately. If 
the technique described above is applied to the second equation, 
say, it appears that the two parameters G and P in the growth 
factor, G exp(in), considered as functions of the wave number, a, 
behave as shown in Figs. 1 and 2. Waves whose length is larger 
than four mesh intervals (a<0.25) are practically neither dif- 
fused nor dispersed. Even waves whose length is comparable to 
two mesh intervals are little diffused, and almost not dispersed 
at all, when the Courant number is close to 1 or 2. This is an 
important feature of the scheme, of which we will take advantage 
when constructing our computational technique. 

A basic FORTRAN code for the scheme described above is 
presented in Fig. 3. 


6 , The role of discontinuities 


Three types of discontinuities must be considered, in order 
of increasing severity: 

1) Gradient discontinuities, across vrtiich all physical 
parameters are continuous but the derivatives of u, a and p jump. 

2) Contact discontinuities, across which s and a junp, 
but p and u are continuous. 
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and 

3) Shocks, across which u, a, s and p jump. 

All discontinuities may be generated at the initial time or 
along boundaries. Typically, gradient discontinuities are the 
head and tail of expansion waves which are produced by a discon- 
tinuity in the acceleration or the velocity of a piston, or by 
the sudden removal of a diaphragm. Contact discontinuities are 
also typically generated at the sudden removal of a diaphragm. 
Shocks are generated at the initial line as compression waves of 
finite strength produced by a discontinuity in the velocity of a 
piston or by the sudden removal of a diaphragm. In addition, 
contact discontinuities are also generated at the collision of 
two shocks of opposite families, or when a shock overtakes anoth- 
er shock of the same family, and shocks are also generated within 
the flow field as a consequence of coalescence of compression 
waves. Finally, gradient discontinuities are generated at the 
intersection of shocks or of shocks and contact discontinuities, 
but in these cases they are sufficiently weak to be neglected. 

Discontinuities are generated and supported by the fact 
that certain signals cannot propagate across them. A computa- 
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tional code which satisfies the same conditions will automatical- 
ly make the discontinuities evident. Otherwise, the discontinui- 
ties will be spread over several mesh points. With an integra- 
tion scheme of second-order accuracy, oscillations will appear on 
either side of the discontinuities, unless special provisions are 
taken to force raonotonicity. 

From a mathematical viewpoint, a function of a real varia- 
ble, continuous with its first derivative (and, if necessary, 
suitably continued over the entire real axis, periodically or 
not) can be considered as the real part of a function, analytic 
on the complex plane (with suitable singularities not lying on 
the real axis). When such functions are smooth (as the physical 
parameters, solutions of Euler's equations, are, when discon- 
tinuities are removed) , their Fourier spectra are dominated by 
very low frequencies. Therefore, if integration schemes of the 
type described in Section 5 are used in such a way that they ap- 
ply only to functions which are continuous with their first 
derivatives, neither diffusion nor dispersion will occur, and the 
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DTrCFL 

CFLrIO. 

LOOP=-1 

13 DO 1 N=1,NC 
DUM=A(N)*GG 
R(1 ,N)=DUM-U(N) 

R(2,N)=DUM+U(N) 

AL(1 ,N)=U(N)-A(N) 
AL(2,N)=U(N)+A(N) 

AL(3,N)=U(N) 

CFL=AMIN 1 ( CFL , STAB*DX/ 

1 (ABS(U(N))+A(N))) 

1 CONTINUE 

DO 2 N=1,NA 
DO 2 L=1,3 

2 RX(L,N)=(R(L.N+1 )-R(L,N))*DDX 
DO 3 N=1,NC 

DO 4 L=1 ,3 
DUM=1 . 

IF(L.EQ.3) DUM=0. 
J=N-(I(L,N)+1)/2 

4 FD(L)=-.5*(AL(L,J)+AL(L,J+1))» 
1 (RX(L,J)-DUM»A(N)*RX(3,J)) 
FD(4)=-A(N)«U(N)»BPRB(N) 
IF(LOOP,EQ.1) GO TO 5 


DO 6 L=1 ,4 
6 F(L,N)=.5»FD(L) 

GO TO 3 
5 DO 8 L=1,3 
J=N-I(L,N) 

8 FN(L,N)=FD(L)-F(L,J) 
FN(4,N)=FD(4)-F(4,N) 

3 CONTINUE 

IF(LOOP.EQ.-I) GO TO 10 
DO 11 N=1,NC 
DO 11 L=1,4 

11 F(L,N)=FN(L,N) 

10 CALL BOUND 

DO 12 N=1,NC 

AT=.5*GD»(F(1,N)+F(2,N)+2. 
1(A(N)»F(3,N)+F(4,N))) 
UT=.5»(F(2,N)-F(1 ,N)) 
ST=F(3,N) 

A(N)=A(N)+AT*DT 

U(N)=U(N)+UT»DT 

12 R(3,N)=R(3,N)+ST*DT 
TIME=TIME+.5*DT 
IF(LOOP.EQ.I) RETURN 
L00P=1 

GO TO 13 


Nomenclature 


GD 

6 

GG 

1/6 

CFL 

At 

DT 

At 

STAB 

Courant number 

R(3,N) 

S 

BPRB 

Ax/A 

N 

n 

NA 

number of intervals 

NC 

NA+1 

AL 

X 

I(L,N) 

sgnCXi^n^ 


Fig. 3 - Basic FORTRAN code for the X-soheme. 



solution will maintain its accuracy and smoothness. 

This philosophy, which the senior author has applied in some 
early papers [13,14,15] is now used again to produce a code which 
aims to be, at the same time, simple, general and robust. The 
same philosophy and code can be extended to multi-dimensional 
problems, whereas codes inspired by [15] are essentially one- 
dimensional. We will show that all that is required to neutral- 
ize the effects of discontinuities and to capture them on a 
single mesh interval are some elementary switches in the x-scheme 
code presented above. 


7. Gradient discontinuities 


Gradient discontinuities are the weakest type of discon- 
tinuities in a flow and, usually, are not considered worthy of a 
special treatment. It is known that gradient discontinuities 
travel along characteristics. If u^ and a^^ are arbitrarily for- 
ced to be discontinuous at (x^.t^), and R 2 jj are discontinuous 
at (x ,t ) but, for t>t^, the discontinuity in R.^ travels along 
the ( u-a)-characteristic issuing from (Xotto^ discon- 
tinuity in travels along the ( u+a)-characteristic Issuing 
from the same point. Therefore, in general, across a gradient 
discontinuity one of the two Riemann variables is always continu- 
ous and differentiable. Its derivative can be approximated using 
a difference between two points on opposite sides of the discon- 
tinuity. The same cannot be said of the other Riemann variable. 
Therefore, indiscriminate application of the x-scheme, as exposed 
in Section 5, in the presence of a gradient discontinuity, can 
produce inaccuracies. Unfortunately, some of these inaccuracies 
can pass undetected, unless one disposes of an exact solution for 
comparison, because they may not affect the local monotonicity of 
the results. It should be kept in mind that, whereas mesh- 
dependent oscillations are a sure index of inaccuracy, monotoni- 
city is not necessarily synonymous with accuracy. 

A gradient discontinuity with an abrupt change in the slope 
of R^, but not of R 2 occurs when u<a. If the discontinuity falls 
somewhere between points A and B, as in Fig. 4a, all x— 
derivatives of R^ and R 2 can be approximated correctly, according 
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to the \-scheme, except at A. If the discontinuity, however, 
IS exactly at A (Fig, 4b), all x-derivatives can be approximated 
correctly. Gradient discontinuities generally move and cross 
over from one interval to another. The calculation benefits from 
the motion, for two reasons: (i) Often the discontinuity is close 
to the position shown in Fig. 4b, and the approximation is almost 
correct; (ii) Also, the error produced at A when the discontinui- 
ty is in the position of Fig, 4a is corrected as soon as A is 
crossed over and starts playing the role of B. 

None of the errors produced by the unconditional use of the 
X-scheme at a gradient discontinuity seems to be catastrophic. 
Indeed, if the discontinuity is the front of an expansion wave, 
the error tends to be diluted, and if it is the front of a 
compression wave, it soon will give rise to a shock wave, which 
has to be treated in a different way. 

Nevertheless, a few disturbing cases of inaccuracies can oc- 
cur which justify a mild intervention on the logic of the 
X-scheme in the vicinity of a gradient discontinuity, as follows: 

1) In the case of Fig. 4a, point A is subject to a pertur- 
bation ccsning from the interior of the advancing wave before the 
wave front reaches A. In turn (and, again, particularly if the 
wave front remains for a while to the right of A) , the error at A 
propagates to the point to the left of it, and the perturbation 
moves at a speed larger than the physical speed of the wave 
front . 


2) The formation of a "precursor" wave in front of the phy- 
sical wave can also be related to dispersion in the numerical 
scheme. As shown above, the x-scheme is non-dispersive for any 
Courant number between 0 and 2 and low wave frequencies. If a 

function has a gradient discontinuity, however, its Fourier spec- 

—2 

trum is of the order of a for increasing o. High order harmon- 
ics have a certain relevance (even as not as strong as in the 
case of jumps, where the spectrian depends on terms of order o”^), 
and we may expect some dispersion to appear. 

3) Stronger effects are expected when the original discon- 
tinuity is stronger; that is, particularly when the perturbation 
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begins as a centered expansion wave. 


The discussion will be resumed in Section 24, where examples 
of unpleasant effects of gradient discontinuities will be given. 
In the present code, we will modify the x-scheme to avoid forma- 
tion of unphysical precursors. To that effect, the derivatives 
of the Riemann variable which has the gradient discontinuity will 
be computed from outside the waves at the points immediately out- 
side; the calculation at points inside will not be modified. A 
local loss in accuracy is inevitable at the points immediately 
outside the waves, but it does not affect the computation at oth- 
er nodes. A similar loss in accuracy is necessary to provide 



(a) Discontinuity between nodes (b) Discontinuity at a node 


Fig. ^4 - Behavior of and Rj at a gradient discontinuity. 


monotonicity in flux difference splitting schemes. 


8. Contact discontinuities 


In Fig. 5, points A and B again bracket the discontinuity. 
Now the velocity u and the pressure p are continuous between A 
and B. Since the entropy is not continuous, it follows that the 
speed of sound a is not continuous either. A simple relation ex- 
ists, however, between the speeds of sound on either side of the 
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contact discontinuity. In the spirit of a discretized code, we 
assume that the same relation exists between the speeds of sound 
at A and B (in other words, we assume that the two points are 
sufficiently close for the pressure to be practically the same at 
both points) : 


a^/ag = exp[6(s^ - Sg)] (25) 

In trying to determine a^ and u^ at A and B, using (13), we see 
that, at A, the (u-a)-wave is always proceeding from the discon- 
tinuity towards its own left, whereas, at b, the (u+a)-wave is 
always proceeding from the discontinuity towards its own right. 
Therefore, f^ is not directly determinable at A and is 

directly determinable at B. From (25), however, we have; 


^tA'^^A “ ^tB^®B "" ~ ®tB^ = ^^^3A “ ^3B^ 

Using the first two equations (13) and (26), and letting u^=Ug, 
f^^ and f2g are easily obtained in the form; 


1A 


1+r 


^2B ■ ^2A ^1B ■* ^1A 


(27) 


where 


Note that the above considerations are correct when both A and B 
are exactly on the contact discontinuity, on opposite sides of 
it, regardless of the flow being subsonic or supersonic. In 
practice, we will force all derivatives at A to be taken from the 
left and all derivatives at B to be taken from the right. Conse- 
quently, will always be computed correctly at A and f^ will be 
computed correctly at B, if the flow at B is subsonic. We must 
allow for some slight inaccuracy in in supersonic flow; al- 
ternatively, we may assume that f^g vanishes identically in su- 
personic flow. 
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9. Shocks 


The strongest discontinuities occur at shocks. Certain 
quantities, however, are conserved across a shock. If the shock 
velocity is denoted by W, then p(u-W), p+p(u-W)^ and (u- 
W)[yp/(y-1)+p(u-W) /2] are the same on both sides of the shock. 
Such properties of conservation (which are one of the many forms 
of the Rankine-Hugoniot conditions) may be useful, once the shock 

p 

velocity, W is known. In an unsteady motion, however, pu, p+pu 
2 

and u[yp/(y- 1 )+Pu /2] are not continuous across a shock. Dif- 
ferentiating them in space across a shock generates high frequen- 
cy waves which can be damped only through the use of artificial 
viscosity, with a consequent loss of accuracy. Other devices 
proposed to improve accuracy, such as the use of adaptive meshes, 
have the intrinsic difficulties of shock fitting and it is diffi- 
cult to extend them to multi-dimensional problems. 


The variables used in the X-scheme are a more natural choice 
for the analysis of a flow which is about to produce a shock. To 
make the point, let us obtain analytical expressions for and 

AR^ across a shock. Here again, we denote by A and B the points 
to the left and to the right of the shock, respectively, and as- 
sume that u>0 and A lies to. the left of the shock, so that the 
high-pressure region (to which B belongs) is to the right of it. 
Let 

M = (-y^)^ (28) 

be the relative Mach number of the shock. Then, 


Au s a 


1-M 
A(1+6)M 


Aa - a 

^ ■ A*^ (1+6)M 


- 1 ] 


(29) 


A simple analysis of (29) for values of M close to 1 shows that 
AR 2 is practically continuous across the shock, whereas Ar.| tends 
to grow much faster with M (the behavior of AR^ and AR 2 func- 
tions of M is illustrated in Figs. 6 and 7). Therefore, the 
X-scheme technique can be used, without having to resort to any 
special provision, until the shock has formed and has taken up 
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some strength. The calculation of at A is correct, not being 
affected by the presence of the shock; it is also acceptable at B 
because of the continuity of across the shock. The calcula- 
tion of R^ at both A and B is correct, because R^ is computed 
from the left at A (where the flow is supersonic) and from the 



right at B (where the flow is subsonic). 

When the shock is about to form, we have indeed all the 
elements which are necessary to detect its formation, and all of 
them are accurate. For the detection of the shock, it is con- 
venient to compute the parameter: 


A (30) 

Between M=1 and M=4, j is a linear function of M to within 1%. 
The exact relation between L and M, obtainable using (29), is 
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( 31 ) 


Z - [(yM^-6)(U6M^)3^^^ 6(M^ - 1) I 

6(U6)M " « 

which can be approximated, for KM<4, by the linear function: 

M = 1 + g r (32) - 

where g is defined by computing (32) at M=4: 

g _ 4S(1 + 6) 

1 16-4+C (16Y -6 )(1+16«)]^^^ (33) 

The detection of the shock is performed by computing Z at every 
point, as defined by (30), and then by applying (32) to obtain an 
approximate value of M. If M is less than 1.1, say, the interval 
is considered as shockless; otherwise, the interval will be la- 
belled as carrying a shock, as it will be explained below. 

When a shock is known to exist, it is computed in a similar 
way, as follows. A first estimate of M is made, using (30) and 
(32). Then a new value, of 2 is obtained from 

s 2 2 - 2' 

where Z’ is computed from (31 )» using the first approximation for 
M. Finally, (32) is applied once more, using Z^, and the evalua- 
tion of M is now accurate to within 10”^. If more accuracy were 
needed, the procedure can be repeated once more. To support our 
conclusions, errors produced by using (32) alone, and by applying 
one or two iterations are plotted in Fig. 8 and labelled 
E^, E2 and E^, respectively. For a given M, the error is the 
difference between the given Z and Z obtained from (31). 


Once the shock Mach number has been obtained , the shock 
velocity, W, follows from (28), and the values at B can be upda- 
ted using the Rankine-Hugoniot conditions, expressed by (29) and 
by 




+ L 


1+0 


V 1+6 /r 


(34) 
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ERROR 


It must be noted that correct values of and are needed at 
B, to be used in the corrector stage for the calculation of the 
point to the right of B. Such values can be obtained by compu- 
ting u^ and s^ numerically, and using the second and third of 
(13). 

Finally, it is very important to note that the calculation 
of the relative Mach number, as outlined above, does not change 
if the high pressure side of the shock lies to the left of the 
discontinuity, provided that and are interchanged and the 
indices A and B are interchanged as well. This is because the 
change in sign in u in the definition of i is compensated by a 
change in sign in M. 

A basic FORTRAN code for the calculation of the shock is 
presented in Fig. 9. 



Fig. 8 - Errors in successive approximations of E(M) using equation (32). 
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L1=3-L 

LL=L-L1 

NN=N+L-1 

NM=NN-LL 

SIGMAr (R (L , NN )-R (L , NM) ) /A( NM) 

TM=1 .+G1»SIGMA 
SM=TM»*2 

SIGMA=2 . »SIGMA-SQRT( (GAMMA»SM-GD) «( 1 .+GD»SM) )+ 
1GD»(SM-1 .))/(GE»TM*GD)+GG 
TM=1 .+G1»SIGMA 
SM=TM**2 

SN=R ( 3 , NM) + ( ALOG( (GAMMA»SM-GD) /GE-GAMMA»ALOG(GE»SM/ 
1(1 .+GD»SM)))/(2.*GD*GAMMA) 

XBT(J)=U(NM)-A(NM)»TM»LL 

R1(L1,NN)=R(L,NN)-2.*(A(NM)*(SM-1.)/(GE*TM)-LL*U(NM)) 

UN=.5«(R(2,NN)-R(1 ,NN)) 

F(L1,NN)=LL»2.*(UN-U(NN))/DT+FL(L,NN) 

F(3,NN)=SN-R(3,NN))/DT 

RETURN 


Nomenclature 


GD 

5 

GG 

1/6 

GE 

1+6 

G1 

g defined by (33) 

SIGMA 

1 

DT 

At 

TM 

M 

R(3,N) 

s 

GAMMA 

Y 

XBT 

*t* 


Fig. 9 - Basie FORTRAN code for shook calculation. 
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1 0 . Shocks in the presence of a contact discontinuity 

The calculation of a shock must be modified when a , contact 
discontinuity belongs to the same interval as the shock. In Fig. 
10, a possible configuration of this type is shown by a plot of 
the entropy distribution. In the same interval, AB, a contact 
discontinuity (C.D.) and a shock coexist. If the technique 
described in Section 9 were applied without modifications, the 
entropy jump across the shock would be automatically assumed as 
the one between A and B, instead of the real one, denoted by (S). 

Let 


= 


‘1A 


- R 


IB 


’B 


(35) 


as it would be evaluated in Section 9, using (30). The proper 
value of 2 should be, instead. 


Z = 



Let 


6 ^ ^s-Sg) 

A = e 


(36) 


(37) 


Then , 


z = 


^o 





(38) 


since the pressure is constant across the contact discontinuity. 

A trial-and-error routine can be used, in which a first 
guess is made of Sg, a is obtained from (37), Z from (38), M from 
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s 



Fig, 10 - Shook and contact discontinuity in the same interval. 


(32), and an error results by comparing the guessed Sg with S, as 
obtained from the right-hand side of (34). 


1 1 . Sonic interval 

Let u > 0 and the flow be subsonic at a point, A and super- 
sonic at the next point, B (Fig. 11). No discontinuity occurs in 
the interval, AB. Nevertheless, the computation must be perfor- 
med with a certain care. Both derivatives of and Rg at A and 
B are computed as prescribed by the x-scheme. What we cannot do 
as usual is the averaging of x.,, as shown in (17), and the 
correction of f^, as shown in (22). Indeed, the prescription 
would carry certain information from A to B and from B to A 
which should not cross over the sonic point. Consequently, 
and the correction to f.| are used with their local values at A 
and B. 
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12. Tracking discontinuities 

We have seen how the basic X“Scheme must be modified locally 
at one end or both ends of an interval containing a discontinuity 
and which additional algebra is needed for certain discontinui- 
ties. We discuss now a simple and general way of tracking any 
number of discontinuities. 

Every point, x = (n-1) Ax (n>^1), is denoted by an index, 
IN(N) (in FORTRAN notation), which encodes the nature of the 
discontinuity or discontinuities, present in the interval ^ the 
right of the point . By a simple device, one single index can 
describe an interval containing more than one discontinuity, 
without ambiguities. Let the discontinuities be denoted by in- 


creasing 

powers 

of 2, IKL) = 2**L (in FORTRAN notation again) 

as follows: 


L=1, 

11=2, 

a shock with the low-pressure region on the 
left. 

L=2. 

11=4, 

a shock with the high-pressure region on the 
left. 

OO 

II 

11=8, 

a contact discontinuity. 
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L=4,5 11=16,32, gradient discontinuities carried by right- 

running characteristics. 

L=6,7 11=64,128, gradient discontinuities carried by left- 

running characteristics. 

In addition, we will use the number 256 (corresponding to L=8) to 
denote an interval in which the flow velocity changes from sub- 
sonic to supersonic. 


In the binary system used by computers, any power of 2 
represented by a 1 , properly located. For example. 


IS 


2 

4 

8 

16 

32 


is represented by 

ti 


10 

100 

1000 

10000 

100000 


etc . 


The coexistence of two or more discontinuities in a single inter- 
val will bring in more than one power of 2. For example, if one 
high-to-low pressure shock,one contact discontinuity and two 
left-running gradient discontinuities exist in the same interval 
(as at the sudden removal of a diaphragm), such an interval will 
be denoted by 

IN(N) =4+8+64+128 

that is, by IN(N)=204, but in binary notation this number will be 

IN(N) = 11001100 


It is clear, thus, that the index still preserves the individual- 
ity of all discontinuities. To reveal each one of them, separa- 
tely, it will be sufficient to perform a logical .AND. between 
the index, IN(N) and the pertinent II(L). In FORTRAN notation, 
if 


IN(N) . AND . IKL) 

is true, the interval contains the discontinuity denoted by 2^. 


To summarize, the new appearance of a discontinuity denoted 
by IKL) with KLj<7, or of a sonic point, denoted by 11(8), en- 
tails a Boolean sum to be coded in FORTRAN as 

IN(N) = IN(N) . OR . IKL) 

The disappearance of a discontinuity or of a sonic point will be 
coded as 

IN(N) = IN(N) . AND . (. NOT . 11(D) 

A discontinuity or a sonic point will be detected by a logical 
IF: 


IF (IN(N) . AND . IKL)) 
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with the expression being true. 

All the operations described above are unambiguous and are 
not time consuming. Note that arithmetic sums cannot be substi- 
tuted for Boolean sums. For example, the coalescence of a shock 
(2) and a contact discontinuity (8) would generate a 10, which 
could also be interpreted as denoting the presence of a shock (2) 
and two other shocks (4+4). 

Note also that two values of L are assigned to either type 
of gradient discontinuity. Such a redundance is necessary only 
to distinguish between the head and the tail of an expansion 
wave, when they belong to the same interval, as it happens when a 
centered expansion wave is generated. Head and tail move at dif- 
ferent speeds; therefore, they must be considered as different 
discontinuities. 

The information concerning discontinuities is stored in two 
arrays, depending on an index, J. The first array contains the 
location, x», of the discontinuity within the interval (0^x*^Ax); 
the second array contains the speed at which the discontinuity 
moves (x^»). Such information is not used in actual updating of 
physical parameters. It is only applied to detect the crossing 
of a discontinuity from one interval to a neighboring one. This 
happens whenever x* becomes negative or larger than ax. The po- 
sition of the discontinuity within the interval is, of course, 
updated using its velocity. 

The velocity of a discontinuity, in turn, is easy to obtain. 
For a gradient discontinuity, the velocity is u+a or u-a; for a 
contact discontinuity, it is u; and the shock velocity, W, is 
computed from (28). 

In the arrays, the information is stored considering each 
interval in order, from N=1 on. As an interval appears which 
contains discontinuities, the pertinent information is stored in 
the J-arrays in order of increasing L. Obviously, when neeeded, 
such information must be retrieved in the same order. 
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13. Treatment of boundaries 


Different problems require different modeling of boundary 
conditions, but the basic requirements at boundaries (regardless 
of their nature) are always the same. If we turn back to Eqs. 
(12), and assume that u>0 and the flow is subsonic, we see that 
at the left boundary f^ and f^ should be computed from outside; 
consequently, they must be considered as unknown. Similarly, at 
the right boundary f^ should be computed from outside and, there- 
fore, must be considered as unknown. If the flow were super son- 
ic, f^, and f^ are all unknown at the left boundary and known 
at the right boundary. 

The modeling of the boundaries provides the codes necessary 
to determine the unknowns without introducing arbitrariness and 
conserving the same computational accuracy as at interior points. 

In the present study, we consider only open (permeable) 
boundaries, as follows: 

1) In strictly one-dimensional problems (flows in ducts 
of constant cross-sectional area) , we assume that both the left 
and the right boundary are reached only by simple waves. In ad- 
dition, we stipulate that the flow outside the left boundary is 
homentropic. Therefore, at the left boundary, f^ = 0 aqd is 
also equal to 0 because is a constant throughout the left- 
running simple wave. For a similar reason, f^ = 0 at the right 
boundary. 

2) In ducts of variable cross-section, we consider only 
nozzles connecting two infinite reservoirs. In the one at the 
left, the flow is at rest, with a uniformly distributed stagna- 
tion temperature and a constant entropy. The second condition 
obviously provides again f^ = 0. To apply the first, let us re- 
call that, in a steady flow, 

a^ = + « u^ (39) 

The time-derivative of (39) at the entrance of the duct must re- 
flect the fact that a„ is constant outside; thus, 

o 
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a * u = 0 (^0) 

at the first computational node. By comparison of (40) with the 
first two Eqs. (13)i we obtain 

(u-a)fi - 2afu 

^ / H ^ \ 


In the reservoir at right, we assume that the pressure is con- 
stant. Such an assumption does not imply that the speed of sound 
is also constant, since the entropy may vary inside the duct, due 
to the formation and motion of shocks. We must now combine (10) 
with the first and the third of (13) to find 

f^ = - f2 - 2f^ + 2aP^/Y (^2) 

3) If the flow is supersonic at the left boundary, we 
assume that it is uniform outside; therefore all f^ equal 0. 


1 4 . Outline of the computational code 

In addition to a main program, which serves as a general or- 
ganizer, the code contains three sets of subroutines. 

Gasdynamical Subroutines 

MARCH - provides the basic integration scheme as described in 
Section 5, with some additional features: 

a) reducing the stepsize if a discontinuity is travel- 
ling over two mesh intervals. 

b) calling TEST to detect sonic points and formation of 

shocks . 

c) calling BOUND to enforce boundary conditions. 

d) calling TRACK to treat discontinuities. 

e) calling XING to provide for shifting of discontinui- 
ties from one interval to the next. 

TEST - detects sonic points and the formation of shocks. 

TRACK - computes contact discontinuities according to Section 8, 
calls SHOCK and contains some instructions for the interaction of 
contact discontinuities and shock and of two wave fronts. 
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SHOCK - computes shocks according to Sections 9 and 10. 

BOUND - computes boundary conditions according to Section 13. 

Logical Subroutines 

XING - provides the logic for the crossing of discontinuities 
from one interval to the next; if necessary, it calls the fol- 
lowing subroutines: 

INSERT - when a new discontinuity appears, 

DELETE - when a discontinuity disappears, 

SWITCH - which accomodates x* and x^* in their general arrays ac- 
cording to the new position of the discontinuity, 

COUNT - which provides the necessary bookkeeping when two discon- 
tinuities from two neighboring intervals are interchanged. 

Initial Condition Subroutines 

XMP, XMQ and XMR - provide different types of initial conditions 
for the examples shown in Sections 15 through 28 and, when 
available, the exact solutions to the problems. 

Additional Logic 

A few additional logical statements are added to the code, 
which are justified by the following considerations. 

If two discontinuities (or more) appear in two consecutive 
intervals, we distinguish two cases: 

1) The discontinuities are the head and the tail of the 
same expansion wave. In this case the calculation may proceed 
according to the logic outlined above. The special provisions 
for the treatment of gradient discontinuities do not conflict 
with each other. 

2) In all other cases, the point bracketed by the 
discontinuities may receive successive Instructions which are in- 
consistent with each other. Instead of going through the cumber- 
some analysis of all possible cases, we prefer to skip the calcu- 
lation of the point. Since all discontinuities move, the point 
will be restored to normality at the next step. When such a si- 
tuation arises, we cut the stepsize in half, in order to reduce 
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the period of neglect. 


Across a contact discontinuity, u does not change. There- 
fore, in principle, the speed of the discontinuity may be taken 
as the value of u at the point to the left of the discontinuity 
or to the right of it. We use the second option in all cases ex- 
cept one, that is, when the contact discontinuity is in the same 
interval as a shock having the high pressure on the left. 

When a contact discontinuity crosses from one interval to 
the next, the entropy at the point which has been crossed over 
must be changed into the value pertinent to the other side of the 
discontinuity and the speed of sound has to be recomputed accor- 
ding to (25). 


15. Test cases. Strictly one-dimensional flows 

To test the efficiency of the logic of the code, we chose 
the following problems in which the mechanics of the flow is not 
complicated by variable area effects: 

1) Collapse of an expansion shock. 

2) Head-on collision of two shocks. 

3) Riemann’s shock tube problem. 

4) Merging of two shocks of the same family. 

5) Decaying shock. 

6) Interaction of a shock and a contact discontinuity. 

7) Interaction of two expansion waves. 

Problems 1,2,3 and 4 have also been treated in [9] (cases 
c,i,f, and i again), and we have tried to use the same data to 
make the comparison between our results and the results of [9] as 
easy as possible. 

Our demonstrative displays consist essentially of two 
figures for each case: a comparison of computed values (denoted 
by O's) and exact values of p (solid lines) at a given time, 
whenever an exact solution is available, and a tracking of 
discontinuities, as computed, in the (x,t)-plane. Discontinui- 
ties are labeled with the same values of L as used in Section 12. 
Occasionally, we will also present sets of lines of constant den- 
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sity (isopycnics) in the (x.t)-plane. 


16. Collapse of a n expan sion shock 

This test, as mentioned in [93, is meant to check the abili- 
ty of the computational code to avoid producing unphysical expan- 
sion shocks. We start assuming that an expansion shock exists, 
and let it collapse. Initially, the flow proceeds from a uniform 
region 1 (defined by a subsonic Mach number, M^, a speed of 
sound, a^ r [y/C 1+6M^) and an entropy, = 0), into a uni- 
form region 2, crossing a shock in reverse. The values in region 
2 are obtained using the Rankine-Hugoniot conditions across the 
expansion shock. The correct flow consists of a left-running ex- 
pansion wave, joining region 1 to a uniform region 4, a contact 
discontinuity, separating region 4 from a uniform region 3, and a 



Fig. 12 - Collapse of an expansion shock. 

very weak shock, separating region 3 from region 2 (Fig. 12). 
Across the wave, a^+«u^ = a^j+su,,, and = P^ + (y/6) ln(a^/a^). 
In addition, u^su^j, P^=Pij=P 2 +ln[ (vM^-6)/( 1+6)] , if M is the Mach 
number of the weak shock. Finally, W being the shock velocity, 
^2 = (u 2 -W)/a 2 and (u 2 ~W)/(u 2 -W) = (1+6M^)/[ (1+6)M^j. These 
equations allow an exact solution to be evaluated by trial-and- 
error. The case presented here corresponds to M^=0.5. Initial- 
ly, the index (IN) for the interval containing the expansion 
shock is set equal to 200 (8+64+128). Fig. 13 shows p(x) at 
t=0.3100. Fig. 14 shows the tracking of the wave fronts (6 and 
7), of the sonic point (8) and of the contact discontinuity (3) 
between t=0 and t=0.3100. Fig. 15 shows the isopycnics for the 
same interval of time. We note: 
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Fig. 73 - Exact and computed density in the collapse of an expansion shock 

(100 intervals). 

1) Perfect adherence of the computed values to the exact 

values. 

2) No overflowing of perturbations across the wave 

front, either in form of oscillations (as in non-upwind schemes 
and in the original x-scheme) or in form of monotonic smearing 
(as in [93) . 

3) No inaccuracies at the sonic point. 

4) Perfect capture of the contact discontinuity. 

5) The weak shock (which, at t=0.3100, is out of the 
picture already) is not tracked, for being too weak. 
Nevertheless, at previous times, a minor jump in p reveals it. 

In runs 102 and 103, the number of intervals has been 

changed to 50 and 25 respectively. Some deterioration, particu- 
larly near the wave front, is expected and it is evident in Figs. 
16 and 17. In Fig. 18, four typical stages in the evolution of 

the flow are presented, as superimposed plots of p(x) for dif- 

ferent values of t, in Run 103. 
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Fig. 15 _ laopycnios in the collapse of an expansion shock. 
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Fig. 18 - Evolution of p(x) in time in collapse of an expansion shock. 
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17. Head-on collision of two shocks 


Here we start assuming the existence of a right-running 
shock, separating a high pressure region (1) of uniform flow from 
a low-pressure region (0) also uniform, and of a left-running 
shock, separating region 0 from another high pressure region, la- 
beled 2, also uniform. The two shocks will collide and continue 
moving in their original directions, at different speeds. A con- 
tact discontinuity will also appear (Fig. 19), For the exact 
solution, the new shock velocities and the physical parameters in 
regions 3 and 4 are obtained by trial-and-error , using the 
Rankine-Hugoniot conditions on the two new shocks and imposing 

P 2 =P 4 , u^su^. The case presented here (Run 201) corresponds to a 
right-running shock defined by M=2 and a left-running shock de- 
fined by M=1.5. The computation is performed on 50 intervals. 
Initially, we set IN=2 in the interval containing the right- 
running shock and IN=4 in the interval containing the left- 
running shock. Fig. 20 shows p(x) at t=0.3107 and Fig. 21 shows 
the tracking of the two shocks. The contact discontinuity is not 
tracked for being too weak; it is, however, detectable as a minor 
step in Fig. 20. Fig. 22 shows some plots of p(x) at successive 
times just before and just after the collision. 


18. Riemann* s problem 

This is a classical case, in which, consequent to an initial 
jump in pressure, a left-running expansion wave, a contact 
discontinuity and a shock are generated (Fig. 23). The same re- 
lations are satisfied as for the previous case of an expansion 
shock, and an exact solution can be obtained accordingly. The 
case presented here corresponds to an initial pressure ratio 
equal to 10 and uniform temperature throughout (so that the en- 
tropy is lower in region 1 than in region 2). Initially, the in- 
terval where the pressure jumps occur has IN=202. In this case 
(Run 301), 100 intervals were used. Fig. 24 shows p(x) at 
t=0.0128. Fig. 25 shows the discontinuities, and Fig. 26 the iso- 
pycnics. The computed shock, in Fig. 24, appears slightly ahead 
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Fig. 20 - Density after the head-on collision of two shocks. 



Fig. 21 - Tracking of shocks in head-on collision. 
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Fig. 22 - Evolution of p(x) in the collision of two shocks. 



Fig. 23 - Riemann's problem. 
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of the exact shock. This is due to an initial uncertainty in the 
location of the shock within the interval and it may or may not 
appear in plots, according to the location of the exact shock 



within the interval at the time being plotted. 


1 9 . Merging of two shocks 

When a shock overtakes another shock of the same family, the 
two merge into a single shock. A contact discontinuity and a 
weak, left-running expansion wave are generated (Fig. 27). The 
exact solution after the merger can be evaluated by the same pro- 
cedure used for the collapse of an expansion shock and for the 
Riemann problem. We started with one shock having M=1.5 and the 
other having M=3. In run 401, we used 100 intervals. The main 
difficulty in this case stems from the narrowness of the expan- 
sion wave which, for a large number of computational steps, 
remains confined to a single interval. In addition, the "left- 
running" wave is actually moving from left to right, so that, 
when the wave front crosses over a mesh point, the point itself 
has to be returned to the unperturbed region. Results of our 
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Fig. 25 - Tracking of discontinuities in Riemann's problem. 
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Fig. 26 - Isopycnics in Riemann's problem. 
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Fig. 27 - Merging of two shocks of the same family. 



Fig. 26 - Density distribution in the merging of two shocks. 
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Fig. 29 - Tracking of discontinuities in merging of tvo shocks. 

calculation are shown in Figs. 28 and 29. It is clear that, 
despite the difficulties just mentioned, the wave moves correctly 
and without smearing. 


Decaying shock 


A fixed shock separates a uniform, supersonic flow from a 
uniform, subsonic flow. A left-running expansion wave propagates 
through the subsonic region and interferes with the shock, which 
loses strength and moves downstream. No perturbation is propaga- 
ted upstream of the shock, since the flow is supersonic in that 
region. Isopycnics are shown in Figs. 30 and 31 to demonstrate 
the effect of mesh size on accuracy. 
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Fig. 31 - Isopyonics for a decaying shock, 50 intervals. 
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21 . Interaction of a shock and a contact discontinuity 


Two different flow patterns may be produced by the interac- 
tion of a shock and a contact discontinuity. If the shock 
penetrates into a region where the speed of sound is lower, the 
shock itself is refracted and another shock emerges, of the oppo- 
site family. A contact discontinuity is also produced. If the 
shock penetrates into a region where the speed of sound is higher, 
the shock is refracted and a centered expansion wave is produced, 
together with a contact discontinuity. The two cases are shown 
in Fig. 32. 


We assume that the flow in regions 1,2, and 3 is given. In 

1 /2 

particular, we assume that U 2 =U 2 = 0 , a^=Y . P 2 =P 3 = 1 f 3^=0; we 
prescribe a^ and compute from 

We also prescribe the Mach number of the shock and obtain all 
pertinent values in region (1) as above. 


In the first case (a 2 >a^), we find the exact solution by 
trial-and-error . The relation to be used, with W^^ = velocity of 
the refracted shock, W^^^= velocity of the reflected shock, are: 


M53 = ”53/^3 


= In [(yM ^2 - «)/(1+6)] 


- 


1/2 

U5 = - l/Mg^) 

CyCyM^^- fiXl + / C( 1 + 6 )M 52 ^ 

s^ = [P^ - Y In (1+6 )m|^y6 


P, = P 5 
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Fig. 32 - Interaction of a shock and a contact discontinuity. 



Fig. 33 - Density after shock-contact discontinuity interaction, second kind 
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?4 = + In [(yM^ 4 - 6)/( 1+6)3 

^14 = - 31^14 

u'4 = [u^(1+6M^4) + W^4(M^4-1)3 / [(1+6)m2^3 

An error results from u^- u*^ and must be made to vanish. 

In the second case a^), a similar procedure is fol- 

lowed. All the equations above are used again, except the last 
three, which must be replaced with: 



1 3 
1 3 
1 3 


X 


Fig. 3^^ - Tracking of discontinuities in shock-contact discontinuity 

interaction, second kind. 
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An error results from 


- a^exp[5(s^-s^)] 
and must be made to vanish. 

In Figs. 33 and 34 we present the results of a calculated 
interaction of the second type. The shock Mach number is 2, the 
contact discontinuity is originally still, and the ratio of the 
speeds of sound between left and right hand side of the discon- 
tinuity is 0.8. Its most interesting feature is the perfect 
description of an extremely narrow expansion wave. 


In Fig. 35 the result of a similar calculation, with a ratio 
of speeds of sound a^/a^sM, is presented. The pattern is now re- 



Fig. 35 - Density after shock-contact discontinuity interaction, first kind. 


versed with respect to Fig. 33 t and a weak shock is reflected. 
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22. Interaction of two expansion waves 


At t=0, two centered expansion waves of opposite families 
are produced at x=0 and x=1. The one at the left is weaker than 
the one at the right. The first is, indeed, all subsonic and all 
contained in the ( x>0)-region. The second has a supersonic part 
which does not appear in the ( x<1 )-region . Heads and tails of 
the waves are shown in Fig. 36, as computed between x=0 and x=1. 
Isopycnics are shown in Fig. 37. They are symmetric up to the in- 
stant of interaction. The subsequent dissymmetry is a consequen- 
ce of the different strengths of the waves. The top of Fig. 37 
shows, on the left, the beginning of a new simple wave. 



X 


Fig* 36 - Tracking of discontinuities in expansion waves Interaction. 


- - 



X 

Fig- 37 - Isopycnios in the interaction of two expansion waves. 



To show the ability of the computational logic to handle a 
large number of discontinuities at the same time, and to take 
proper care of their interactions, we have considered the flow 
produced by the simultaneous breaking of two diaphragms. Initi- 
ally, the flow appears as the succession of two similar patterns, 
of the type shown in Fig. 25. An interaction of the type shown 
in Fig. 31 is the first to occur, complicated by the motion of 
the flow behind the moving shock. A shock-contact discontinuity 
interaction will eventually occur but is not shown in the figure 
(Fig. 38). Isopycnics are shown in Fig. 39. The shock to the 
right is too weak to appear, since the number of level lines is 
automatically controlled by the program to avoid excessive 
clustering of lines. 
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Fig. 39 


Isopycnios for two simultaneous Riemann problems. 
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24 . Numerical exe rcis es o n gradient d is continuities 

In Section 7i we have discussed the numerical problems rela- 
ted to the presence of gradient discontinuities, in general 
terms. The code described in Section 14 can be easily modified 
to delete the tracking of gradient discontinuities, in order to 
support our arguments with some quantitative experiments. 

We begin considering a centered expansion wave originated at 
x=0.1682, t=-0.2804, by a piston moving backwards with a velocity 
equal to -0.6. The exact solution is computed at t=0 (note that 
in this case both a and u are linear within the wave, at a con- 
stant t) . The computation is then started at t=0, using the 
X-scheme as originally described in Section 5, at all points. 
The result is almost perfect (Fig. 40); the example does not sug- 
gest any need for special trackings or modifications to the 
scheme. 

The same computation is repeated (Fig. 41) reducing the com- 
putational stepsize by a factor of 10. A precursor disturbance 
appears, although it is not large enough to be considered as a 
threat to the accuracy of the computation. Note that in Figs. 40 
and 41 the speed of sound, not the density, is plotted; the 
results, however, are the same, qualitatively. This exercise is 
made to support our arguments in Section 5. It must be noted 
that, in a complicated flow, the stepsize may have to be drasti- 
cally reduced at certain points, and all other pxjints are affec- 
ted; therefore, the possibility of dispersion of wave fronts due 
to low Courant numbers is not a remote one. 

More drastic effects are shown in Figs. 42 and 43 (where 
density again is on the ordinate). Here, the expansion wave is 
centered at t=0, x=0 and the motion of the piston is so fast that 
the flow becomes sonic at the left boundary. The calculation 
performed with a tracking of the gradient discontinuity produces 
results which fall exactly on the solid line. Without tracking, 
the results show a sizeable precursor and a consequent degenera- 
tion over the entire wave (Fig. 42). The effect is, of course, 
enhanced when the stepsize is reduced by a factor of 10 (Fig. 
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43 ) . 


The next figure (Fig. 44) shows the result of the calcula- 
tion of the same Riemann problem presented in Fig. 24. The 
deterioration of the expansion wave is evident and does not re- 
quire any comment. 

A similar deterioration appears in Fig. 45, which should be 
compared with Fig. 28. At the time considered in the figure, the 
expansion wave is barely encompassing four intervals. Without 



Fig. MO - Speed of sound in a simple expansion wave, no tracking. 


tracking, the wave is spread over 16 intervals. 


25 . More test cases. Quasi-one-dimensional flows 

The following problems were considered as tests of the accu- 
racy of the numerical analysis, rather than of the efficiency of 
the logic. To make the flow as far from uniform as possible, we 
consider ducts of variable cross-sectional area. The area, A(x), 
is defined by imposing a Mach number distribution between x=0 and 
x=1 and determining A(x) as satisfying the area-Mach number rela- 
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X 

Fig. ^5 - Density for merging shocks without tracking of gradient discontinuity. 


tion for a steady, homentropic flow: 

A = (1 + (i»3) 

In turn, the Mach number distribution is either linear (prescri- 
bed by the values of M at x=0 and x=1) or parabolic (prescribed 
by equal values of M at x=0, x =1 and by its value at x=1/2). 

If we assume that the flow starts from a given initial state 
(for example, from rest, with a sudden drop in pressure at x=1), 
but that it remains homentropic, and that the duct proceeds from 
an infinite cavity containing a gas at rest, the flow inside the 
duct will approach a steady state asymptotically, whose Mach num- 
ber distribution should be the same as the one initially chosen 
to define A(x). We have, thus, an asymptotic exact solution with 
which the computed flow may be compared. 

If the original Mach number crosses from subsonic to super- 
sonic, and the prescribed exit pressure is properly chosen, the 
asymptotic steady state will no longer be homentropic. A steady 
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shock will appear in the divergent part of the duct, followed by 
a subsonic region. Once more, we have at our disposal the exact 
solution, which is easily obtained by prescribing the shock loca- 
tion and calculating the jump in the critical area. From the 
shock to the exit, M can be computed from (43), altering A by a 
constant factor. It will also be easy to find which exit pres- 
sure must be prescribed, consistent with a given shock location. 


26. Subsonic , homentropic flows 

1) If M is prescribed to go from 0.2 to 0.8 linearly (con- 
verging duct) , an asymptotic steady state is reached very ra- 
pidly. If the distribution is the opposite one (diverging duct), 
it takes much longer to reach a steady state. These well-known 
facts are confirmed by the isopycnics of the two flows (Figs. 46 



and 47). The obvious steady state results are not worth showing. 

2) If M is prescribed as a parabolic curve, with entry and 
exit values equal to 0.3 and a maximum value equal to 0.8, the 
isopycnic plot appears as in Fig. 48 and the asymptotic values 
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X 

Fig. *(8 - Isopyonics for parabolic duct. 
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3 



Fig. 51 - Area distributions in parabolic ducts. 


(denoted by O’s) fall very well on top of the exact steady state 
solution (solid line, in Fig. 49). 

Difficulties appear if the lower value of M is 0.2 instead 
of 0.3 and the upper value is 0.9 instead of 0.8. When starting 
from rest, after an initial phase in which the computed values 
become very close to the steady state solution, a second phase 
follows with some drifting in the upper subsonic region; eventu- 
ally, the computed flow develops a small supersonic region, fol- 
lowed by a weak shock (Fig. 50). The failure of the code, which, 
at a first sight, could be attributed to lack of accuracy in the 
transonic region, is to be found in a lack of accuracy in the en- 
try portion of the duct. The curve of A(x) in this case appears 
as shown in Fig. 51; for comparison, A(x) is also shown for the 
preceding case, 0.3 < M < 0.8, and for the case, 0.3 < M < 0.9. 
It seems that a quasi-one-diraensional approach is, in the present 
instance, inadequate to approximate the real behavior of a flow 
in a duct whose cross-sectional area varies so rapidly. Numeri- 
cal errors in the calculation of the entry point may produce a 
slight error in stagnation enthalpy and a consequent change in 
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TIME 


the critical area. In addition, the cross-sectional area remains 
flat and very close to 1 over a large portion of the duct; there- 
fore, the change in the critical area may be sufficient to repla- 
ce the prescribed duct with a duct having a throat exactly at 
x=0.5. The Mach number distribution shown in Fig. 50 is indeed 
consistent with such a geometry. Any other Mach number distribu- 
tion generating channels without such a steep area distribution 
can be computed with accuracy. It is interesting to note that 
all these computations with parabolic Mach number distributions 
are extremely long, because the initial perturbation is very 
small and it requires a large number of waves, travelling back 
and forth along the channel , in order to produce a fast flow in 
the middle of it. 



Fig. 52 - Isopycnios in Laval nozzle with a shock. 


27. Transonic flows with a shock 

We present in Figs. 52, 53 and 54 the results obtained by 
prescribing a duct geometry consistent with an isentropic linear 
Mach number distribution between M=0.5 and M=2.0, and an exit 
pressure consistent with the presence of a shock at x=0.63. As 
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in the preceding cases, the flow is generated by the sudden drop 
of the exit pressure to the final value, the flow being original- 
ly at rest. The isopycnic plot clearly reveals the original re- 
ceding expansion wave and its reflection as a compression wave, 
which eventually coalesces into a shock. 

Finally, Fig. 55 shows the isopycnic pattern in the same 
duct, assuming that, for t<0, the flow in the duct was steady, 
with a shock somewhere between x=0.63 and x=1. At t=0 the exit 
pressure is suddenly lifted to the value of the previous case. 
The shock is pushed backwards to the new stable position, x=0.63. 
The final Mach number distribution is not shown, being the same 
as in Fig. 54. 


28 . An exercise on contact discontinuities 

To conclude our numerical experiments, we present calcula- 
tions made on two flows in a parabolic nozzle. In both cases, 
the flows start from rest. The initial entropy, however, is not 
constant. It is prescribed to be equal to 0 to the left of the 

throat and to have another value, s^ = (ln6)/6 to the right of 

^ 1/2 1/2 

it. The speed of sound equals and 8y at the left and 

right of the throat, respectively. If the speed of sound at the 
exit is dropped to a fraction, a, of its original value, a flow 
is generated and the contact discontinuity eventually moves 
throughout the divergent portion of the duct and exits from it. 
At a certain stage of the evolution, according to the value of 6, 
the contact discontinuity may separate a subsonic flow from a su- 
personic flow, or vice versa. The case of “=0.7, ^=1.2 is shown 
in Figs. 56, 57 and 58, the case of “=0.7, ^=0.7 in Figs. 59, 60 
and 61 . A history of the evolution of the Mach number is dep- 
icted in Figs. 56 and 59, to show how the transonic transition 
occurs. The purpose of these calculations is to show that no 
symptoms of errors appear, regardless of the nature of the tran- 
sition . 
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Fig. 57 - Isopyonics in nozzle with contact discontinuity, a^>a^. 



X 

Fig. 58 - Tracking of discontinuities in nozzle, a 2 >a.|. 
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Fig. 60 - Isopyonios in nozzle with contact discontinuity, 
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Fig. 61 - Tracking of discontinuities in nozzle, a^<ay 


Concluding remarks 

An improved method for the calculation of one-dimensional 
flows has been presented. It combines a simple and efficient 
version of the A -scheme with tracking of discontinuities. The 
latter is needed to identify points where minor departures from 
the basic integration scheme are applied to prevent infiltration 
of numerical errors. Such a tracking is obtained via a systema- 
tic application of Boolean algebra. It is, therefore, very effi- 
cient. Fifteen examples of flow computations have been presented 
and discussed in detail. The results are exceptionally good. 
All discontinuities are captured within one mesh interval. Ef- 
forts are currently being made to extend the method to multidi- 
mensional flows. 
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